COO - single cell RNA sequencing data analysis

Bioinformatics in Neuroscience course, 2022

The course used the CoCalc server: https://cocalc7.science.uu.nl/

Some good practice
To keep everything ordered, make some folders for saving the data and the figures. You will receive a warning if you have already made these folders

dir.create("write/")
dir.create("figures/")

1.1 Load libraries

Load the necessary libraries to run our pipeline based on the well-known Seurat pipeline. You can find detailed information on Seurat as well as tutorials on its website (https://satijalab.org/seurat/) or on Github (https://github.com/satijalab/seurat)

We are going to use Seurat v3, which is well established. Note that there is a recent beta release (v4) with additional functionalities

# Load the libraries
library(Seurat,verbose = FALSE)
library(ggplot2,verbose = FALSE)
library(RColorBrewer,verbose = FALSE)
library(useful,verbose = FALSE)
library(readr,verbose = FALSE)
library(hdf5r,verbose = FALSE)

# check if the following are on
library(dplyr,verbose = FALSE)
library(stringr,verbose = FALSE)

1.2 Load the data

We will use the 10x Genomics data generated from the adult human brain. Since 4M cells are too many, we are using a ~5k subset. The data is stored in two seperate files

# This contains the actual data as a cell x gene matrix. Each cell in this matrix contain information on how many reads we have for a given gene in a given cell
count_data <- read_rds("../rawdata/A46_10X176_8_count_data.rds")

# the dimensions for a sanity check. The first dimension shows the number of features (in this case, genes) while the second indicates teh number of cells
dim(count_data)

# This contains the metadata. Each column contains information for each cell in teh dataset.
metadata <-  read_rds('../rawdata/A46_10X176_8_metadata.rds')

# what do we have in the metadata?
colnames(metadata)

1.3 Create Seurat object

Here, we will create a "seurat object" which has a different structure than a data table. Seurat object contains the data as a matrix in a separate slot, but can also have 'features' that are linked to genes or cells. For instance cell types for each cell can be stored in a separate slot. In addition, we can store additional information e.g. results of differential expression test or coordinates for each cell after a given dimensionality reduction (e.g. pca)

'min.cells' is to filter genes: only genes detected in min of these number of cells will be kept 'min.features' is to filter cells: only cells that have a minimum total of given number of features will be kept Thus, these two provide initial quality control of the data

dataset <- CreateSeuratObject(counts = count_data, min.cells = 1, min.features = 1000, meta.data = metadata,
                          project = "Human_brain")

visualise the number of genes and features, which are already provided in the dataset

head(dataset@meta.data)

The first visual quality control is to look at the number of cells (nCount_RNA) and number of genes (nFeature_RNA) per cells. These features are calculated while generating the Seurat object. For instance, the number of RNA molecules per cell is stored as 'nCount_RNA', which is located in dataset@meta.data$nCount_RNA

1.4 Visualise basic metrics

plot1 <- VlnPlot(object = dataset, features = c("nCount_RNA", "nFeature_RNA"),pt.size = .01)
plot2 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 | plot2

2.1 Standard pre-processing workflow

set a name that will be used later for the files you save. Anything goes. You can add your own name

name = "my_analysis_name"

2.1.1 Calculate additional QC metrics

mRNA coding for the Ribosomal subunit proteins is abundant (not to be confused with rRNA). Usually, a high ribosomal RNA percentage indicates production of a lot of proteins, and is very high in dividing cells or some secretory cells that need to constantly produce proteins. However, if most of the mRNA (>30-50%) that we detect is ribosomal, it means that we wont learn much from this cell, and that we should exclude it from the analysis.

# Start by generating  QC metrics additional to the no of genes/features
dataset <- PercentageFeatureSet(dataset,pattern='^MT-', col.name = "percent.mt") #Select all features starting with MT-, which are the Mitochondrial genes
dataset <- PercentageFeatureSet(dataset,pattern='RP(S|L)', col.name = "percent.ribo") # Select all that starts with RPS or RPL

Plot and inspect. You can remove the ots in case they cause too much moise

2.1.2 Visualise QC metrics

options(repr.plot.width=12, repr.plot.height=6)
VlnPlot(object = dataset, features = c("percent.mt", "percent.ribo"))

You can save them all as a pdf file to use in presentations, reports etc

pdf(paste0(name,"_Vlnplot_stats_all.pdf"),height = 6,width = 12)
VlnPlot(object = dataset, features = c("nCount_RNA", "nFeature_RNA","percent.mt", "percent.ribo"))
# If there are too many points for your liking, you can turn them off as follows:
#VlnPlot(object = dataset, features = c("nCount_RNA", "nFeature_RNA","percent.mt", "percent.ribo"),pt.size = 0, cols = "blue")
dev.off()

Visualise multiple metrics for each cell as a measure of the number of reads (in other words; UMI, or unique RNA molecules). This helps use to evaluate the quality of the data in multiple ways

plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.ribo")
plot1 | plot2

2.2 Filter

Decide on the filtering parameters that you want to use. This is not always objective! However, you can see the 'outliers' in the data rather well

#Now filter
dataset <- subset(x = dataset, subset = percent.mt < 5 & nCount_RNA < 90000 & percent.ribo < 2)

dataset

Now get rid of the mitochondrial genes which we do not need for downstream analysis.

dataset <- subset(x = dataset[-grep(pattern = "^MT-", x = rownames(dataset), value = FALSE)])
dataset 

3.1 Normalisation

Standard flow of Seurat is replaced by a single command. However, it is good to see this part to learn each step

## normalise the data
dataset <- NormalizeData(object = dataset, normalization.method = "LogNormalize", scale.factor = 10000)

The log transformed RNA counts from each cell is now normalised using a size factor. In other words, the sum of counts for each cells is more or less the same. We can now compare different cells to each other

3.2 Detection of variable genes across the single cells

## Here, we select the top 1,000 highly variable genes (Hvg) for downstream analysis.
dataset <- FindVariableFeatures(object = dataset, selection.method = 'mean.var.plot', mean.cutoff = c(0.0125, 3), dispersion.cutoff = c(0.5, Inf))
length(x = VariableFeatures(object = dataset)) #3084
## Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(dataset), 10)

## Plot
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(dataset)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
#plot1
plot2

# Generate a vector with the top 1000 hvgens
hv.genes <- head(rownames(HVFInfo(object = dataset)), 1000)

Scale the data

Scale the values for each features. The result is in units of 'standard deviation' where the mean expression is zero. If a gene is expressed 1 standard deviation more than the mean across all cells, it will have a value of 1. During this process, we can also 'regress out' technical effects. For instance, the number of features and the mitochondrial percentage.

dataset <- ScaleData(object = dataset, features = hv.genes, vars.to.regress = c("nFeature_RNA") ) 

4.1 Perform dimensionality reduction

You can run PCA, t-SNE and UMAP in simple one-line codes... The parameters are NOT random, but we will ignore this for the sake of time

# Dimention reduction
dataset <- RunPCA(object = dataset, features = hv.genes, verbose = FALSE,npcs = 50)

dataset <- RunUMAP(object = dataset, reduction = "pca", dims = 1:50, verbose = FALSE)
dataset <- RunTSNE(object = dataset, dims = 1:50, verbose = FALSE)

4.2 Examine and visualize PCA results a few different ways

4.2.1 Dimplots

First, let's have a look at the first three dimensions. For this, we will take advantage of some of hte metadata generated by the authors: Supercluster information

options(repr.plot.width=12, repr.plot.height=6)
DimPlot(dataset,cols = 'Spectral',group.by = 'Supercluster',pt.size = 1,reduction = 'pca',dims = c(1,2))
DimPlot(dataset,cols = 'Spectral',group.by = 'Supercluster',pt.size = 1,reduction = 'pca',dims = c(3,4))

4.2.2 List of genes that influence PCs that most in both (+/-) directions

You see the top genes that contribute to the respective principle components

print(dataset[["pca"]], dims = 1:5, nfeatures = 5)

4.2.3 Visualise genes that influence PCs the most

VizDimLoadings(dataset, dims = 1:2, reduction = "pca")

4.2.3 Heatmaps to visualise genes that influence PCs

options(repr.plot.width=12, repr.plot.height=9)
DimHeatmap(dataset, dims = 1:4, cells = 500, balanced = TRUE, ncol = 2)

4.3 Determine statistically significant principal components

How many principle components do we need? There are some practical guidelines that will help with your analysis

4.3.1 JackStrawPlot()

We can do this with JackStrawPlot() function of Seurat. This is powerful, but time consuming. Running it altogether slows the server considerably. Also, we cannot visualise all different PCs (only up to 25). We will skip it for now (you can run it during the practice if you want to) and instead look at PCs using heatmaps

dataset <- JackStraw(object = dataset, num.replicate = 100)
dataset <- ScoreJackStraw(object = dataset, dims = 1:20)
JackStrawPlot(object = dataset, dims = 1:20)

4.3.2 Heatmaps

We will plot the top genes from selected PCs and see how they differe between cells.

options(repr.plot.width=18, repr.plot.height=9)
DimHeatmap(dataset, dims = c(1,2,5,10,
                             15,20,25,30,
                             35,40,45,50), cells = 250, balanced = TRUE, ncol = 4)

4.2.2 Elbowplot as an alternative

# Plots the standard deviations (or approximate singular values if running PCAFast) of the principle components for easy identification of an elbow in the graph. This elbow often corresponds well with the significant dims and speeds up running than Jackstraw
# In this example, it looks like the elbow would fall around PC 11, then 25 and 35 again. After this, there is little change. 
options(repr.plot.width=12, repr.plot.height=6)
ElbowPlot(object = dataset, reduction = "pca",ndims = 50)

5.1 Find clusters

here, we will use graph-based clustering to identify groups of cells that are similar, which we will call 'clusters'

5.1.1 Use embeded 'louvain' algorithm for clustering

#SNN Graph Construction
dataset <- FindNeighbors(object = dataset, dims = 1:35) #cell the pcs based on Jackstraw and Elbow plots
#Find clusters using the louvain algorithm
dataset <- FindClusters(object = dataset, resolution = 1) # changing the resolution will change the size/number of clusters! c(0.6, 0.8, 1)

5.1.2 Visualize QCs

Let's have a look at the QC stats agian. As default, this time the information will be displayed per cluster

# Statistics
options(repr.plot.width=12, repr.plot.height=6)
VlnPlot(object = dataset, features = c("nCount_RNA", "nFeature_RNA"),pt.size = 0)

5.2 Cell type identification

5.2.1 Visualise lineage specific genes curated from the literature

We provide some markers from the literature as a reference for identifying cell types. You can find more of these in Siletti et al 2022 Bioarxiv manuscript

ELAVL2, TUBB3: Neurons PLP1: Oligodendrocytes GFAP: Astrocytes ITGAM: Microglia PECAM1: Endothelial cells

5.2.1.1 t-SNE plots for visualising markers at single cell level

options(repr.plot.width=18, repr.plot.height=12)
ColorRamp <- colorRampPalette(rev(brewer.pal(n = 7,name = "RdYlBu")))(100) # This is how to generate a color key
FeaturePlot(object = dataset, features = c("ELAVL2","TUBB3","PLP1","GFAP","ITGAM","PECAM1"),reduction = 'tsne',cols = ColorRamp, ncol =3)

Neurons are very diverse! We can use additonal markers to classify them functionally as 'Glutamatergic' and 'GABAergic'

SLC17A5, SLC17A6: Glutamatergic neurons GAD1, GAD2, SLC32A1: GABAergic neurons

options(repr.plot.width=18, repr.plot.height=12)
ColorRamp <- colorRampPalette(rev(brewer.pal(n = 7,name = "RdYlBu")))(100) # This is how to generate a color key
FeaturePlot(object = dataset, features = c("ELAVL2","GAD1","GAD2","SLC17A5","SLC17A6","SLC32A1"),reduction = 'tsne',cols = ColorRamp, ncol =3)

5.2.1.2 Violin plot for cluster level average

options(repr.plot.width=18, repr.plot.height=6)
VlnPlot(object = dataset, features = c("ELAVL2","TUBB3","PLP1","GFAP","ITGAM","PECAM1"),pt.size = 0, ncol =2) #, slot = 'counts'
# you can plot raw counts as well

options(repr.plot.width=18, repr.plot.height=6)
VlnPlot(object = dataset, features = c("ELAVL2","GAD1","GAD2","SLC17A5","SLC17A6","SLC32A1"),pt.size = 0, ncol =2) #, slot = 'counts'
# you can plot raw counts as well

You can use marker gene information to name each cluster as

cell classes e.g. oligodendrocytes, neurons or marker-defined cell types e.g. GlutamatergicNeuron_0

5.2.2 Compare our clusters with the metadata

Your clusters are saved in the metadata, as part of the seurat object, with the name "seurat_clusters". We compare this to the metadata from the authors, which will help you decide which clusters belong to which cell type

plot21 <- DimPlot(object = dataset, reduction = 'tsne',group.by = "seurat_clusters")
plot22 <- DimPlot(object = dataset, reduction = 'tsne',group.by = "Supercluster")
plot23 <- DimPlot(object = dataset, reduction = 'tsne',group.by = "Clusters")

#plot, adjusting the size of the graph
options(repr.plot.width=18, repr.plot.height=6)
plot21 | plot22
options(repr.plot.width=12, repr.plot.height=6)
plot23

We can visualise the number of cells labeled with each label using the table() function

table(dataset@meta.data[['Supercluster']])

5.2.3 Visualise how different labels compare to each other

To compare them, we will do a trick. First, generate a dataframe with teh information you need, then use the table() function of make a confusion matrix

df <- as.data.frame(cbind(dataset@meta.data[['Supercluster']],dataset@meta.data[['seurat_clusters']]))
rownames(df) <- colnames(dataset)
colnames(df) <- c('Supercluster','seurat_clusters') #Supercluster Clusters

head(table(df)[,1:10]) #head() automatically looks at teh first 6 rows. We also shorthen the columns to 10

Quick hetmap to visualise which cell types belong to which cluster

heatmap(table(df),xlab ='seurat_clusters')

6.1 find marker genes for each cluster

6.1.1 Start with a single cluster

We often want to see what makes a cluster special, or different than the rest. You can achieve this here

cluster1.markers <- FindMarkers(object = dataset, ident.1 = 0, min.pct = 0.25) #, ident.2 = c(0, 3) to compare
head(x = cluster1.markers, n = 5)

6.1.2 Find markers for every cluster compared to all remaining cells, report only the positive ones

This is extremely useful to define what these cells are! However, you will see that a 'global comparison' is not always as accurate as you would expect

library(dplyr)
dataset.markers <- FindAllMarkers(object = dataset, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
# we have just calculated the marker genes for each cell type

We can select the top genes for each cluster. Here, we will select the top 2 genes

dataset.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)  %>% head(n=10)
# Takes a bit of time

save the data as separate sheets in an excel file

library(openxlsx)
## Create Workbook object and add worksheets
wb <- createWorkbook()
## Add the cluster info
for ( n in unique(dataset.markers$cluster) ) {
  addWorksheet(wb, n)
  writeData(wb, n, dataset.markers[dataset.markers$cluster == n,], startCol = 1, startRow = 1, rowNames = TRUE, colNames = TRUE)
}
#save
saveWorkbook(wb, paste(name,"_cell_clust_diff_genes.xlsx",sep=""), overwrite = TRUE)

6.2 DoHeatmap generates an expression heatmap for given cells and features.

In this case, we are plotting the top 20 markers (or all markers if less than 20) for each cluster.

dataset.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) -> top10

options(repr.plot.width=12, repr.plot.height=12)
#dataset <- ScaleData(object = dataset, features = rownames(top10$gene)) #,vars.to.regress = "percent.mt"
DoHeatmap(object = dataset, features = top10$gene[1:100],slot =  'scale.data') + NoLegend()

More information: https://satijalab.org/seurat/

https://satijalab.org/seurat/v3.1/visualization_vignette.html

Extra code for the practical

subsetting neuronal clusters '0' and '1'

small <- subset(dataset, idents = c('0','1'))

small

DimPlot(small, reduction = 'tsne')

You can proceed from step 3.2 to repeat the analysis of a subset of hte dataset that contains clusters 1 and 2

THE END